# -*- coding: utf-8 -*-

from pylab import *
from scipy import signal

def reversed_bits(x, width):
    """Генерирует последовательность битов числа x в обратном порядке.
    
    В отличие от binary_repr, биты представляются целыми числами, 
    а не строками.
    """
    
    for bit in reversed(binary_repr(x, width)):
        yield bit

def bin_to_dec(bin_repr):
    """Переводит целое неотрицательное число из двоичной в десятичную систему 
    счисления.
    
    bin_repr: строка двоичного представления числа.
    """
    
    dec = 0
    num_width = len(bin_repr)
    for i in xrange(num_width):
        if bin_repr[i] == '1':
            dec |= 1 << (num_width - i - 1)
    return dec

def seq_to_had(i, width):
    """Переводит номер функции Уолша при упорядочении по частости (Уолшу) 
    в её номер при упорядочении по Адамару."""
    
    gray_i = i ^ (i >> 1)
    had_binary_repr = list(reversed_bits(gray_i, width))
    had = bin_to_dec(had_binary_repr)
    return had

def dya_to_had(i, width):
    """Переводит номер функции Уолша при упорядочении по Пели (диадическое)
    в её номер при упорядочении по Адамару."""
    
    had = list(reversed_bits(i, width))
    return bin_to_dec(had)

def fwht(x, ordering = 'sequency'):
    """Выполняет быстрое преобразование Уолша c заданной сортировкой."""
    
    # вычисление спектра
    tx = do_fwht(x)
    # сортировка спектра согласно указанному порядку
    otx = sort_spectr(tx, ordering) * 1.0 / tx.size
    return otx

def ifwht(x, ordering = 'sequency'):
    """Выполняет обратное быстрое преобразование Уолша с заданной 
    сортировкой."""
    
    tx = do_fwht(x)
    return sort_spectr(tx, ordering)

def do_fwht(x):
    """Выполняет быстрое преобразование Уолша с сортировкой функций по Адамару.
    
    Эта функция выполняет само преобразование. Функции fwht лишь сортируют 
    результат работы этой функции.
    """
    
    n = int(log2(x.size))
    tx = array(x)
    # выполняем преобразование
    #for i in xrange(1, n + 1):
    for i in xrange(1, n + 1):
        numOfBlocks = 1 << (n - i)
        for b in xrange(numOfBlocks):
            halfLenOfBlock = 1 << (i - 1)
            for k in xrange(halfLenOfBlock):
                # индексы элементов для выполнения бабочки
                id1 = b * (1 << i) + k
                id2 = id1 + halfLenOfBlock
                # выполняем бабочку
                tmp1, tmp2 = tx[id1], tx[id2]
                tx[id1] = tmp1 + tmp2
                tx[id2] = tmp1 - tmp2
    return tx

def sort_spectr(tx, ordering):
    """Сортирует спектр tx, отсортированный по Адамару.
    
    tx: результат работы функции do_fwht
    """
    
    n = int(log2(tx.size))
    # выбираем нужную функцию преобразования индексов
    if ordering == 'sequency':
        to_had = lambda i: seq_to_had(i, n)
    elif ordering == 'hadamard':
        # спектр изначально отсортирован по Адамару
        return array(tx)
    elif ordering == 'dyadic':
        to_had = lambda i: dya_to_had(i, n)
    else:
        raise ValuerError

    # сортируем
    otx = array(tx)
    for i in xrange(tx.size):
        otx[i] = tx[to_had(i)]
    return otx
